# include "scene.h"
# include <stdio.h>
# include <math.h>

extern char **parse_str(const char *str, char c, unsigned short start, unsigned short end, unsigned short *nsubs);
extern double str2num(const char *str, short nsubs);

static char r_conf_grid(struct grid_config *gc, struct grid *g);
static char th_conf_grid(struct grid_config *gc, struct grid *g);
static char z_conf_grid(struct grid_config *gc, struct grid *g);

# ifndef READ_GRID_FILES
static int nearest_int(double x);
# endif

extern void conf_grid(struct grid_config *gc, struct grid *g, char *ecode)
{
 ecode[0] = ecode[1] = ecode[2] = 'o';

 ecode[0] = r_conf_grid(gc, g);
 ecode[1] = th_conf_grid(gc, g);
 ecode[2] = z_conf_grid(gc, g);

 return;
}



static char r_conf_grid(struct grid_config *gc, struct grid *g)
{
 unsigned short izone, isec, i, n;
 FILE *ptr;
 char *str, **s, fname[MAX_FNAME_LEN];
 double temp;

#ifdef READ_GRID_FILES
/* Opening the r-grid file to load the grids from. */
 sprintf(fname, "%02d.%s.dat", RUN_CNT-1, R_FILE);
 fprintf(g->lptr, "Reading grids from the file '%s'.\n", fname);

 ptr = (FILE *)fopen(fname, "r");

/* Determining the number of char in first-line. */
 for(n = 0; getc(ptr) != '\n'; n++);
 str = (char *)malloc(sizeof(char)*(n+1));
 rewind(ptr);

/* Reading first line of the file to determine number of grids and other grid parameters for various sections. */
 for(i = 0; i < n; i++) str[i] = getc(ptr);
 str[n] = '\0';
 s = parse_str(str, ' ', 1, Nsec+3, &n);
 free(str);

 if(n != Nsec+3)
     {
      fprintf(g->lptr, "ERROR:: In 'r_conf_grid()': the r-grid file '%s' seems to be currupt.\n", fname);
      return 'r';
     }

/* Assigning number of grids for various sections. */
 for(isec = 0; isec < Nsec; isec++)
     {
      g->nr[isec] = (unsigned short)str2num(s[isec], 1);
      free(s[isec]);
     }

 g->nr_nozzle = (unsigned short)str2num(s[Nsec], 1); free(s[Nsec]);
 g->nr_gap = (unsigned short)str2num(s[Nsec+1], 1); free(s[Nsec+1]);
 if((n = (unsigned short)str2num(s[Nsec+2], 1)) != NrG)
     {
      fprintf(g->lptr, "ERROR:: In 'r_conf_grid()': In the file '%s', number of ghost cells = %d != NrG = %d.\n", fname, n, NrG);
      return 'r';
     }

 free(s[Nsec+2]); free(s); fclose(ptr);

#else

# ifdef r_UNIF_GRID
 fprintf(g->lptr, "ERROR:: In 'r_conf_grid()': No rule defined for grid generaion in r-direction.\n");
 return 'r';
# else

 fprintf(g->lptr, "Non-uniform grids are being generated in r-direction.\n");

/* The Algorithm used to generate the r-grids is as follows:

 1. The top-section-1's grids are generated first and then these are copied to section-0 and 2.
 2. The top-section-1 is divided into 6 zones (0 ~ 5):
	* zone-0: 0 ~ r_NOZZLE
	* zone-1: r_NOZZLE ~ r0_CYL
	* zone-2: r0_CYL ~ r1_MID1
	* zone-3: r1_MID1 ~ r1_MID2
	* zone-4: r1_MID2 ~ r1_OUTLET
	* zone-5: r1_OUTLET ~ r1_CYL
 3. "zone-1 & zone-2" and "zone-4 & zone-5" are mirror images.
 4. zones 0,1, and 5 are created independently at first.
 5. Then, zones 2 and 4 are created as a mirror image of zones 1 and 5, respectively.
 6. Then, zone-3 is created to connect the zone-2 and zone-4 smoothly. In this process, r1_OUTLET and r1_CYL changes. */

 gc[0].nzone_r = N0rZONE;
 gc[1].nzone_r = N1rZONE;
 gc[2].nzone_r = N2rZONE;

 gc[1].r[0] = (double)0;
 gc[1].r[1] = r_NOZZLE;
 gc[1].r[2] = r0_CYL;
 gc[1].dr[0] = dr_CENTER;
 gc[1].dr[1] = dr_NOZZLE; /* Modified. */
 gc[1].dr[2] = dr0_CYL; /* Modified. */

 gc[1].r[5] = r1_OUTLET; /* Modified. */
 gc[1].r[6] = r1_CYL; /* Modified. */
 gc[1].dr[5] = dr1_OUTLET;
 gc[1].dr[6] = dr1_CYL; /* Modified. */

 for(g->nr[1] = 0, izone = 0; izone <= 5; izone++)
    {
     gc[1].nr[izone] = nearest_int(2*(gc[1].r[izone+1]-gc[1].r[izone])/(gc[1].dr[izone]+gc[1].dr[izone+1]));

     if( gc[1].nr[izone] < 0)
          {
           fprintf(g->lptr, "ERROR:: In 'r_conf_grid()': gc[1].nz[%d] = %d < 0.\n", izone, gc[1].nz[izone]);
           return 'r';
          }
     else if(gc[1].nr[izone] > 0)
          {
           gc[1].dr[izone+1] = 2*(gc[1].r[izone+1]-gc[1].r[izone])/gc[1].nr[izone] - gc[1].dr[izone];
           gc[1].r_cd[izone] = (gc[1].dr[izone+1]-gc[1].dr[izone])/(gc[1].nr[izone]-1);
          }
     g->nr[1] += gc[1].nr[izone];

     if(izone == 1) izone = 4;
    }


 gc[1].r_cd[2] = -gc[1].r_cd[1];
 gc[1].r[3] = r1_MID1; /* Modified. */
 gc[1].nr[2] = nearest_int( (gc[1].dr[2]/gc[1].r_cd[2]-0.5)\
              *( sqrt(1+8*(gc[1].r[3]-gc[1].r[2])*gc[1].r_cd[2]/((2*gc[1].dr[2]-gc[1].r_cd[2])*(2*gc[1].dr[2]-gc[1].r_cd[2]))) - 1 ) );
 if(gc[1].nr[2] < 0)
     {
      fprintf(g->lptr, "ERROR:: In 'r_conf_grid()': gc[1].nr[2] = %d < 0.\n", gc[1].nr[2]);
      return 'r';
     }
 gc[1].dr[3] = gc[1].dr[2] + (gc[1].nr[2]-1)*gc[1].r_cd[2];
 gc[1].r[3] = gc[1].r[2] + 0.5*gc[1].nr[2]*(gc[1].dr[2]+gc[1].dr[3]);
 g->nr[1] += gc[1].nr[2];

 gc[1].r_cd[4] = -gc[1].r_cd[5];
 gc[1].r[4] = r1_MID2; /* Modified. */
 gc[1].nr[4] = nearest_int( (gc[1].dr[5]/gc[1].r_cd[4]+0.5)\
               *( 1 - sqrt(1 - 8*(gc[1].r[5]-gc[1].r[4])*gc[1].r_cd[4]/((2*gc[1].dr[5]+gc[1].r_cd[4])*(2*gc[1].dr[5]+gc[1].r_cd[4]))) ) );
 if(gc[1].nr[4] < 0)
     {
      fprintf(g->lptr, "ERROR:: In 'r_conf_grid()': gc[1].nr[2] = %d < 0.\n", gc[1].nr[2]);
      return 'r';
     }
 gc[1].dr[4] = gc[1].dr[5] - (gc[1].nr[4]-1)*gc[1].r_cd[4];
 gc[1].r[4] = gc[1].r[5] - 0.5*gc[1].nr[4]*(gc[1].dr[4]+gc[1].dr[5]);
 g->nr[1] += gc[1].nr[4];

 temp = gc[1].r[4];
 gc[1].nr[3] = nearest_int(2*(gc[1].r[4]-gc[1].r[3])/(gc[1].dr[4]+gc[1].dr[3]));
 gc[1].r[4] = gc[1].r[3] + 0.5*gc[1].nr[3]*(gc[1].dr[4]+gc[1].dr[3]);
 gc[1].r_cd[3] = (gc[1].dr[4]-gc[1].dr[3])/(gc[1].nr[3]-1);
 g->nr[1] += gc[1].nr[3];
 temp = gc[1].r[4] - temp; /* Calculating the shift in r[4]. */
 gc[1].r[5] += temp; /* Shifting the outlet radius. */
 gc[1].r[6] += temp; /* Shifting the radius of the outer cylinder. */

 fprintf(g->lptr, "Radii of the OUTLET and the main-cylinder, both, are increased by: %E mm\n", temp*1000);

 g->nr_nozzle = gc[1].nr[0];
 g->nr_gap = gc[1].nr[5];

/* Copying the grid paramaters from Section-1 to Section-0. */
 gc[0].r[0] = (double)0;
 gc[0].r[1] = gc[1].r[1];
 gc[0].r[2] = gc[1].r[2];
 gc[0].dr[0] = gc[1].dr[0];
 gc[0].dr[1] = gc[1].dr[1];
 gc[0].dr[2] = gc[1].dr[2];
 gc[0].nr[0] = gc[1].nr[0];
 gc[0].nr[1] = gc[1].nr[1];
 gc[0].r_cd[0] = gc[1].r_cd[0];
 gc[0].r_cd[1] = gc[1].r_cd[1];
 for(g->nr[0] = 0, izone = 0; izone < gc[0].nzone_r; g->nr[0] += gc[0].nr[izone++]);

/* Copying the grid paramaters from Section-1 to Section-2. */
 gc[2].r[0] = (double)0;
 gc[2].r[1] = gc[1].r[6]-gc[1].r[5];
 gc[2].dr[0] = gc[1].dr[5];
 gc[2].dr[1] = gc[1].dr[6];
 gc[2].nr[0] = gc[1].nr[5];
 gc[2].r_cd[0] = gc[1].r_cd[5];
 for(g->nr[2] = 0, izone = 0; izone < gc[2].nzone_r; g->nr[2] += gc[2].nr[izone++]);

# endif

#endif

 return 'o';
}


static char th_conf_grid(struct grid_config *gc, struct grid *g)
{
 unsigned short i, n;
 FILE *ptr;
 char **s, *str, fname[MAX_FNAME_LEN+1];

#ifdef READ_GRID_FILES
/* Opening the theta-grid file to load the grids from. */
 sprintf(fname, "%02d.%s.dat", RUN_CNT-1, TH_FILE);
 fprintf(g->lptr, "Reading grids from the file '%s'.\n", fname);
 ptr = (FILE *)fopen(fname, "r");

/* Determining the number of char in first-line. */
 for(n = 0; getc(ptr) != '\n'; n++);
 str = (char *)malloc(sizeof(char)*(n+1));
 rewind(ptr);

/* Reading first line of the file to determine number of grids and other grid parameters for various sections. */
 for(i = 0; i < n; i++) str[i] = getc(ptr);
 str[n] = '\0';
 s = parse_str(str, ' ', 1, 2, &n);
 free(str);

 if(n != 2)
     {
      fprintf(g->lptr, "ERROR:: In 'th_conf_grid()': the theta-grid file '%s' seems to be currupt.\n", fname);
      return 't';
     }

/* Assigning number of grids for various sections. */
 g->nth = (unsigned short)str2num(s[0], 1);
 free(s[0]);

 if((n = (unsigned short)str2num(s[1], 1)) != NthG)
     {
      fprintf(g->lptr, "ERROR:: In 'th_conf_grid()': In the file '%s', number of ghost cells = %d != NthG = %d.\n", fname, n, NthG);
      return 't';
     }

 free(s[1]); free(s); fclose(ptr);

#else

# ifdef th_UNIF_GRID
 fprintf(g->lptr, "Uniform grids are being generated in theta-direction.\n");

 g->dth = dTHETA;
#ifdef nTH_COMP
 g->nth = nTH_COMP;
#else
 g->nth = (int)((double)360/g->dth); //Total number of cells in theta direction.
 if(g->nth%2 != 0) fprintf(g->lptr, "No. of cells in theta-direction must be even!!\nIncreasing it by one to nth = %d\n", ++(g->nth));

 g->dth = 360.0/(double)g->nth; //Re-defining "d_theta": in degrees.
#endif
 g->dth *= acos(0.0)/90.0; //"d_theta" in radians.

 if(g->nth < 2*NthG+1) // No. of CELLs chosen for computation keeping in mind symmetry of the problem.
     {
      fprintf(g->lptr, "ERROR:: in 'th_conf_grid()': Defined number of computational cells = %d. However, it should be more than '%d'. Redefine the parameter 'nTH_COMP' to anything more than '%d', re-compile, and rum the program again.\n", g->nth, 2*NthG, 2*NthG);
      return 't';
     }

# else

 fprintf(g->lptr, "ERROR:: in 'th_conf_grid()': No rule defined to generate non-uniform grids in theta-direction.\n");
 return 't';
# endif

#endif

 return 'o';
}

 //[Modify both the conf_grid.c and grid_gen.c accordingly]

static char z_conf_grid(struct grid_config *gc, struct grid *g)
{
 short izone;
 unsigned short n, i, isec;
 FILE *ptr;
 char *str, **s, fname[MAX_FNAME_LEN+1];
 double dist0, dist1, dist;

#ifdef READ_GRID_FILES
/* Opening the z-grid file to load the grids from. */
 sprintf(fname, "%02d.%s.dat", RUN_CNT-1, Z_FILE);
 fprintf(g->lptr, "Reading grids from the file '%s'.\n", fname);
 ptr = (FILE *)fopen(fname, "r");

/* Determining the number of char in first-line. */
 for(n = 0; getc(ptr) != '\n'; n++);
 str = (char *)malloc(sizeof(char)*(n+1));
 rewind(ptr);

/* Reading first line of the file to determine number of grids and other grid parameters for various sections. */
 for(i = 0; i < n; i++) str[i] = getc(ptr);
 str[n] = '\0';
 s = parse_str(str, ' ', 1, Nsec+1, &n);
 free(str);

 if(n != Nsec+1)
     {
      fprintf(g->lptr, "ERROR:: In 'z_conf_grid()': the z-grid file '%s' seems to be currupt.\n", fname);
      return 'z';
     }

/* Assigning number of grids for various sections. */
 for(isec = 0; isec < Nsec; isec++)
     {
      g->nz[isec] = (unsigned short)str2num(s[isec], 1);
      free(s[isec]);
     }

 if((n = (unsigned short)str2num(s[Nsec], 1)) != NzG)
     {
      fprintf(g->lptr, "ERROR:: In 'z_conf_grid()': In the file '%s', number of ghost cells = %d != NzG = %d.\n", fname, n, NzG);
      return 'z';
     }

 free(s[Nsec]); free(s); fclose(ptr);

#else

# ifdef z_UNIF_GRID
  fprintf(g->lptr, "Uniform grids are being generated in z-direction.\n");
  flag = -1;
# else

 fprintf(g->lptr, "Non-uniform grids are being generated in z-direction.\n");

/* For the top section-1 */
 gc[1].nzone_z = N1zZONE;

 gc[1].h[0] = (double)0; gc[1].dh[0] = dh1_BOT;
 gc[1].h[1] = h1_MID1; gc[1].dh[1] = dh1_MID1;
 gc[1].h[2] = h1_MID2; gc[1].dh[2] = dh1_MID2;
 gc[1].h[3] = h1_CYL;  gc[1].dh[3] = dh1_TOP;

 g->nz[1] = 0;

 for(izone = gc[1].nzone_z-1; izone >= 0; izone--) /* Starting from top to bottom */
     {
      gc[1].nz[izone] = nearest_int(2*(gc[1].h[izone+1]-gc[1].h[izone])/(gc[1].dh[izone+1]+gc[1].dh[izone]));

      if(gc[1].nz[izone]  < 0)
          {
           fprintf(g->lptr, "ERROR:: In 'z_conf_grid()': gc[1].nz[%d] = %d < 0.\n", izone, gc[1].nz[izone]);
           return 'z';
          }
      else if(gc[1].nz[izone] > 0)
          {
           gc[1].dh[izone] = 2*(gc[1].h[izone+1]-gc[1].h[izone])/gc[1].nz[izone] - gc[1].dh[izone+1];
           gc[1].h_cd[izone] = (gc[1].dh[izone+1]-gc[1].dh[izone])/(gc[1].nz[izone]-1);
          }

      g->nz[1] += gc[1].nz[izone];
     }


/* For the bottom SECTION-0 */
 gc[0].nzone_z = N0zZONE; /* The top-zone is the mirrot image of the bottom-zone of the SECTION-2. */

 gc[0].h[3] = h0_CYL; gc[0].dh[3] = gc[1].dh[0]; gc[0].h_cd[2] = -gc[1].h_cd[0];
 dist = h0_CYL - h0_MID2;
 dist0 = (double)0;
 for(i = 0; ; i++)
        {
         dist1 = 0.5*(i+1)*(2*gc[0].dh[3] - i*gc[0].h_cd[2]);
         if((dist-dist0)*(dist-dist1) <= 0)
            {
             if(0.5*(dist0+dist1) > dist)
                 {
                  gc[0].h[2] = h0_CYL - dist0;
                  gc[0].nz[2] = i;
                 }
             else{
                  gc[0].h[2] = h0_CYL - dist1;
                  gc[0].nz[2] = i+1;
                 }
             break;
            }
          dist0 = dist1;
         }

 if(gc[0].nz[2] < NzG) /* At least NzG number of cells should be mirror. It is required for the application of BC. */
         {
          fprintf(g->lptr, "ERROR:: In 'z_conf_grid()': gc[0].nz[2] = %d < NzG = %d.\n", gc[0].nz[2], NzG);
          return 'z';
         }

 gc[0].dh[2] = gc[0].dh[3] - (gc[0].nz[2]-1)*gc[0].h_cd[2];

 gc[0].h[0] = (double)0; gc[0].dh[0] = dh0_BOT;
 gc[0].h[1] = h0_MID1; gc[0].dh[1] = dh0_MID1;

 g->nz[0] = gc[0].nz[2];
 for(izone = gc[0].nzone_z-2; izone >= 0; izone--) /* Starting from top to down */
     {
      gc[0].nz[izone] = nearest_int(2*(gc[0].h[izone+1]-gc[0].h[izone])/(gc[0].dh[izone+1]+gc[0].dh[izone]));

      if(gc[0].nz[izone]  < 0)
          {
           fprintf(g->lptr, "ERROR:: In 'z_conf_grid()': gc[0].nz[%d] = %d < 0.\n", izone, gc[0].nz[izone]);
           return 'z';
          }
      else if(gc[0].nz[izone] > 0)
          {
           gc[0].dh[izone] = 2*(gc[0].h[izone+1]-gc[0].h[izone])/gc[0].nz[izone] - gc[0].dh[izone+1];
           gc[0].h_cd[izone] = (gc[0].dh[izone+1]-gc[0].dh[izone])/(gc[0].nz[izone]-1);
          }
      g->nz[0] += gc[0].nz[izone];
     }


/* For the extended SECTION-2 */
 gc[2].nzone_z = N2zZONE; /* The top-zone is the mirror image of the bottom-zone of the SECTION-2. */

 gc[2].h[4] = h2_CYL; gc[2].dh[4] = gc[1].dh[0]; gc[2].h_cd[3] = -gc[1].h_cd[0];
 dist = h2_CYL - h2_MID2;
 dist0 = (double)0;
 for(i = 0; i < 50; i++)
        {
         dist1 = 0.5*(i+1)*(2*gc[2].dh[4] - i*gc[2].h_cd[3]);
         if((dist-dist0)*(dist-dist1) <= 0)
            {
             if(0.5*(dist0+dist1) > dist)
                 {
                  gc[2].h[3] = h2_CYL - dist0;
                  gc[2].nz[3] = i;
                 }
             else{
                  gc[2].h[3] = h2_CYL - dist1;
                  gc[2].nz[3] = i+1;
                 }
             break;
            }
          dist0 = dist1;
         }

 if(gc[2].nz[3] < NzG) /* At least NzG number of cells should be mirror. It is required for the application of BC. */
         {
          fprintf(g->lptr, "ERROR:: In 'z_conf_grid()': gc[2].nz[3] = %d < NzG = %d.\n", gc[2].nz[3], NzG);
          return 'z';
         }

 gc[2].dh[3] = gc[2].dh[4] - (gc[2].nz[3]-1)*gc[2].h_cd[3];

 gc[2].h[0] = (double)0; gc[2].dh[0] = dh2_UNIF_BOT;  gc[2].h_cd[0] = (double)0; gc[2].nz[0] = nz2_UNIF_BOT;
 gc[2].h[1] = nz2_UNIF_BOT*dh2_UNIF_BOT; gc[2].dh[1] = dh2_UNIF_BOT;

/* These two will be modified. */
 gc[2].h[2] = h2_MID1; gc[2].dh[2] = dh2_MID1;

 gc[2].nz[1] = nearest_int(2*(gc[2].h[2]-gc[2].h[1])/(gc[2].dh[1]+gc[2].dh[2]));
 gc[2].nz[2] = nearest_int(2*(gc[2].h[3]-gc[2].h[2])/(gc[2].dh[2]+gc[2].dh[3]));

 gc[2].dh[2] = (2*(gc[2].h[3]-gc[2].h[1]) - (gc[2].nz[1]*gc[2].dh[1]+gc[2].nz[2]*gc[2].dh[3]))/(gc[2].nz[1]+gc[2].nz[2]);
 gc[2].h[2] = gc[2].h[3] - 0.5*gc[2].nz[2]*(gc[2].dh[2]+gc[2].dh[3]);

 gc[2].h_cd[1] = (gc[2].dh[2]-gc[2].dh[1])/(gc[2].nz[1]-1);
 gc[2].h_cd[2] = (gc[2].dh[3]-gc[2].dh[2])/(gc[2].nz[2]-1);

 for(g->nz[2] = 0, izone = 0; izone < gc[2].nzone_z; g->nz[2] += gc[2].nz[izone++]);

# endif

#endif

 return 'o';
}

# ifndef READ_GRID_FILES
static int nearest_int(double x)
   {
    short flag = 1;

    if(x < (double)0)
        {
         x = -x;
         flag = -1;
        }

    return flag*((x-(int)x < (int)x+1-x) ? (int)x : (int)x+1);
   }
# endif
